Drivers of Dissimilarity Walk Through (Bottom temp, balanced Bray Curtis Dissimilarity)

Here, we explore bottom temperature and fishing as potential drivers of spatial beta diversity across trawl regions.

library(data.table)
library(ggplot2)
library(lme4)
library(lmerTest)
library(MuMIn)

###Pull in datatable with dissimilarities, reg characteristics, fishing, and temperature

#from local (will have to change)
dissimilarities_temp_fishing_regstats <- readRDS(here::here("output","distance_decay","dissimilarities_temp_fishing_regstats.rds"))

###Palette for Plotting Palette for plotting all 37 survey units (Prep for eventual plots)

survey_unit.list <- levels(dissimilarities_temp_fishing_regstats[,factor(survey_unit)])

palette_37 <- c(
  "#5A5156", #AI
  "#F6222E", #CHL
  "#F8A19F", #DFO-NF
  "#16FF32", #DFO-QCS
  "#DF00DB", #BITS-1
  "#DB8EDA", #BITS-4
  "#325A9B", #EBS
  "#3283FE", #EVHOE
  "#FEAF16", #FR-CGFS
  "#1C8356", #GMEX-Summer
  "#C4451C", #GOA
  "#85660D", #GRL-DE
  "#B0009F", #GSL-N
  "#BF79B8", #GSL-S
  "#1CBE4F", #ICE-GFS
  "#782AB6", #IE-IGFS
  "#90AD1C", #MEDITS
  "#6B003A", #NAM
  "#A75B00", #NEUS-Fall
  "#E3B072", #NEUS-Spring
  "#02E8B6", #NIGFS-1
  "#97E7D5", #NIGFS-4
  "#B00068", #Nor-BTS-3
  "#00B9E3", #NS-IBTS-1
  "#95E2F4", #NS-IBTS-3
  "#B3CE73", #NZ-CHAT
  "#689500", #NZ-ECSI
  "#AAF400", #NZ-WCSI
  "#AA0DFE", #PT-IBTS
  "#FA0087", #S-GEORG
  "#DEA0FD", #SCS-Summer
  "#FCEF88", #SEUS-fall
  "#A59405", #SEUS-spring
  "#FCE100", #SEUS-summer
  "#C075A6", #WCANN
  "#BDCDFF", #ZAF-ATL
  "#003EFF"  #ZAF-IND
)

color_link <- data.table(survey_unit = survey_unit.list,hex = palette_37)

Add names for plotting

name_helper <- data.table(Survey_Name_Season = c("Aleutian Islands",
                                    "Baltic Sea Q1",
                                    "Baltic Sea Q4",
                                    "Chile",
                                    "Newfoundland",
                                    "Queen Charlotte Sound",
                                    "Eastern Bering Sea",
                                    "Bay of Biscay",
                                    "English Channel",
                                    "Gulf of Mexico",
                                    "Gulf of Alaska",
                                    "Greenland",
                                    "N Gulf of St. Lawrence",
                                    "S Gulf of St. Lawrence",
                                    "Iceland",
                                    "Irish Sea",
                                    "Mediterranean",
                                    "Namibia",
                                    "NE US Fall",
                                    "NE US Spring",
                                    "N Ireland Q1",
                                    "N Ireland Q4",
                                    "Norway",
                                    "N Sea Q1",
                                    "N Sea Q3",
                                    "Chatham Rise",
                                    "E Coast S Island NZ",
                                    "W Coast S Island NZ",
                                    "Portugal",
                                    "S Georgia Straight",
                                  "Scotian Shelf",
                                  "SE US Fall",
                                  "SE US Spring",
                                  "SE US Summer",
                                  "W Coast US",
                                  "Atlantic Ocean ZA",
                                  "Indian Ocean ZA"),
                          survey_unit = c(
                                  "AI",        
                                  "BITS-1",    
                                  "BITS-4",    
                                  "CHL",       
                                  "DFO-NF",    
                                  "DFO-QCS",   
                                  "EBS",       
                                  "EVHOE",     
                                  "FR-CGFS",   
                                  "GMEX-Summer",
                                  "GOA",       
                                  "GRL-DE",    
                                  "GSL-N",     
                                  "GSL-S",     
                                  "ICE-GFS",   
                                  "IE-IGFS",   
                                  "MEDITS",    
                                  "NAM",       
                                  "NEUS-Fall", 
                                  "NEUS-Spring",
                                  "NIGFS-1",   
                                  "NIGFS-4",   
                                  "Nor-BTS-3", 
                                  "NS-IBTS-1", 
                                  "NS-IBTS-3", 
                                  "NZ-CHAT",   
                                  "NZ-ECSI",   
                                  "NZ-WCSI",   
                                  "PT-IBTS",   
                                  "S-GEORG",   
                                  "SCS-SUMMER",
                                  "SEUS-fall", 
                                  "SEUS-spring",
                                  "SEUS-summer",
                                  "WCANN",     
                                  "ZAF-ATL",   
                                  "ZAF-IND"   
                          ))

color_link <- color_link[name_helper, on = "survey_unit"]

###First, we built mixed models with raw temperature predictors as fixed effects and survey as random slope and intercept - Mean bottom temp 12 months before survey - Max bottom temp 12 years before survey - Min bottom temp 12 years before survey - Bottom temp seasonality 12 months before survey - Heterogeneity (SD) of mean bottom temp 12 months before survey - Heterogeneity (SD) of max bottom temp 12 months before survey - Heterogeneity (SD) of min bottom temp 12 months before survey - Heterogeneity (SD) of bottom temp seasonalty 12 months before survey

#from local (will have to change)
balanced_dissimilarity_sbt_model_results <- readRDS(here::here("output","balanced_dissimilarity_sbt_model_results.rds"))

setorder(balanced_dissimilarity_sbt_model_results,-AICc)

balanced_dissimilarity_sbt_model_results[Scaled == F,]

The standard deviation of mean temperature and the overall mean temperature perform equally well. However, neither explain very much variation at all (R^ of 0.7 and 0.2)

The SD of mean tempereature performs best, but explains little variation. *Note that this plot is just to help visualize data and patterns, it does not plot predicted values from the LME but rather just simple linear models of dissimilarity ~ temp.

ggplot(data = dissimilarities_temp_fishing_regstats) +
  labs(x = "SD of mean bottom temperature",  y = "Bray Curtis balanced dissimilarity") +
  geom_point(aes(x = yearly_mean_bypoint_SD, y = bray_curtis_dissimilarity_balanced_mean, color = survey_unit), alpha = 0.3) +
  geom_line(aes(x = yearly_mean_bypoint_SD, y = bray_curtis_dissimilarity_balanced_mean, color = survey_unit), stat="smooth", method = "lm") +
  scale_color_manual(values = c(palette_37)) +
  geom_smooth(aes(x = yearly_mean_bypoint_SD, y = bray_curtis_dissimilarity_balanced_mean), method = "lm", se = F, color  = "black") +
  theme_classic()

#facet helpful for visualization
ggplot(data = dissimilarities_temp_fishing_regstats) +
  labs(x = "SD of mean bottom temperature",  y = "Bray Curtis balanced dissimilarity") +
  geom_point(aes(x = yearly_mean_bypoint_SD, y = bray_curtis_dissimilarity_balanced_mean, color = survey_unit), alpha = 0.3) +
  geom_line(aes(x = yearly_mean_bypoint_SD, y = bray_curtis_dissimilarity_balanced_mean, color = survey_unit), stat="smooth", method = "lm") +
  scale_color_manual(values = c(palette_37)) +
  geom_smooth(aes(x = yearly_mean_bypoint_SD, y = bray_curtis_dissimilarity_balanced_mean), method = "lm", se = F, color  = "black") +
  facet_wrap(~survey_unit, scales = "free") +
  theme_classic() + theme(legend.position = "null")

The mean temperature performs equally well, but explains little variation. *Note that this plot is just to help visualize data and patterns, it does not plot predicted values from the LME but rather just simple linear models of dissimilarity ~ temp.

ggplot(data = dissimilarities_temp_fishing_regstats) +
  labs(x = "Mean bottom temperature (˚C)",  y = "Bray Curtis balanced dissimilarity") +
  geom_point(aes(x = yearly_mean_bypoint_avg, y = bray_curtis_dissimilarity_balanced_mean, color = survey_unit), alpha = 0.3) +
  geom_line(aes(x = yearly_mean_bypoint_avg, y = bray_curtis_dissimilarity_balanced_mean, color = survey_unit), stat="smooth", method = "lm") +
  scale_color_manual(values = c(palette_37)) +
  geom_smooth(aes(x = yearly_mean_bypoint_avg, y = bray_curtis_dissimilarity_balanced_mean), method = "lm", se = F, color  = "black") +
  theme_classic()

#faceted helpful for visualization
ggplot(data = dissimilarities_temp_fishing_regstats) +
  labs(x = "Mean bottom temperature (˚C)",  y = "Bray Curtis balanced dissimilarity") +
  geom_point(aes(x = yearly_mean_bypoint_avg, y = bray_curtis_dissimilarity_balanced_mean, color = survey_unit), alpha = 0.3) +
  scale_color_manual(values = c(palette_37)) +
  facet_wrap(~survey_unit, scales = "free") +
  geom_smooth(aes(x = yearly_mean_bypoint_avg, y = bray_curtis_dissimilarity_balanced_mean), method = "lm", se = F, color  = "black") +
  theme_classic() +
  theme(legend.position = "null")

###Then, we looked at relative temperature predictors (scaled within a survey region, so, is this year warm or cool for this region) - Relative mean bottom temp 12 months before survey - Relative max bottom temp 12 years before survey - Relative min bottom temp 12 years before survey - Relative bottom temp seasonality 12 months before survey - Relative heterogeneity (SD) of mean bottom temp 12 months before survey - Relative heterogeneity (SD) of max bottom temp 12 months before survey - Relative heterogeneity (SD) of min bottom temp 12 months before survey - Relative heterogeneity (SD) of bottom temp seasonalty 12 months before survey

balanced_dissimilarity_sbt_model_results[Scaled == T,]

The relative mean tempereature performs best, but explains 0 variation. *Note that this plot is just to help visualize data and patterns, it does not plot predicted values from the LME but rather just simple linear models of dissimilarity ~ temp.

ggplot(data = dissimilarities_temp_fishing_regstats) +
  labs(x = "Mean bottom temperature scaled",  y = "Bray Curtis balanced dissimilarity") +
  geom_point(aes(x = yearly_mean_bypoint_avg.s, y = bray_curtis_dissimilarity_balanced_mean, color = survey_unit), alpha = 0.3) +
  geom_line(aes(x = yearly_mean_bypoint_avg.s, y = bray_curtis_dissimilarity_balanced_mean, color = survey_unit), stat="smooth", method = "lm") +
  scale_color_manual(values = c(palette_37)) +
  geom_smooth(aes(x = yearly_mean_bypoint_avg.s, y = bray_curtis_dissimilarity_balanced_mean), method = "lm", se = F, color  = "black") +
  theme_classic()

#faceted helpful for visualization
ggplot(data = dissimilarities_temp_fishing_regstats) +
  labs(x = "Mean bottom temperature scaled",  y = "Bray Curtis balanced dissimilarity") +
  geom_point(aes(x = yearly_mean_bypoint_avg.s, y = bray_curtis_dissimilarity_balanced_mean, color = survey_unit), alpha = 0.3) +
  geom_line(aes(x = yearly_mean_bypoint_avg.s, y = bray_curtis_dissimilarity_balanced_mean, color = survey_unit), stat="smooth", method = "lm") +
  scale_color_manual(values = c(palette_37)) +
  geom_smooth(aes(x = yearly_mean_bypoint_avg.s, y = bray_curtis_dissimilarity_balanced_mean), method = "lm", se = F, color  = "black") +
  facet_wrap(~survey_unit, scales = "free") +
  theme_classic() + theme(legend.position = "null")

###Fishing Pressure Alone

#model ranking for fishing only models

#from local (will have to change)
balanced_dissimilarity_fishing_model_results <-readRDS(here::here("output","balanced_dissimilarity_fishing_model_results.Rds"))

setorder(balanced_dissimilarity_fishing_model_results, -AICc)

balanced_dissimilarity_fishing_model_results

Best model only includes survey unit (probably because of differences in baseline dissimilarity across regions)

###Now we’ll look at both fishing and bottom temperature as potential predictors

#model rankings (temp and fishing)
#from local (will have to change)
balanced_dissimilarity_sbt_fishing_model_results <- readRDS(here::here("output", "balanced_dissimilarity_sbt_fishing_model_results.Rds"))

setorder(balanced_dissimilarity_sbt_fishing_model_results, -AICc)
balanced_dissimilarity_sbt_fishing_model_results

Any model with survey as a fixed effect performs better than models without survey as a fixed effect. Also, there is an interaction between fishing and survey in 4 of the top 5 performing models.

Significant positive interaction (higher dissimilarity at higher fishing) - Newfoundland

Significant negative interaction (higher dissimilarity at lower fishing) - Eastern Bering Sea - France - Gulf of Alaska - Gulf of St. Lawrence S - Northeast US Fall - Ireland 4th quarter - North Sea (1st and 3rd quarter) - West Coast and East Coast of South Island of New Zealand - South Georgia - Scotian Shelf Summer Survey

Plot just dissimilarity ~ relative fishing pressure coefficients of each region (Like figure 1b) *Note, may be better to include survey as a fixed effect, but below we include as random slope and intercept instead.

#lmer with fishing as fixed and survey as random slope and intercept
#fishing and random slope and intercept for survey
allreg_dissimilarity_fishing_mean_mod <- lmer(bray_curtis_dissimilarity_balanced_mean ~ summed_tonnes_scaled_byreg + (1 + summed_tonnes_scaled_byreg|survey_unit), data = dissimilarities_temp_fishing_regstats)
                 

 # see group coefficients and confidence intervals
fishing_model_coefs_reduced <- data.table(transform(as.data.frame(ranef(allreg_dissimilarity_fishing_mean_mod)), lwr = condval - 1.96*condsd, upr = condval + 1.96*condsd))
#https://stackoverflow.com/questions/69805532/extract-the-confidence-intervals-of-lmer-random-effects-plotted-with-dotplotra


#ONLY SLOPES
fishing_model_coefs_reduced <- fishing_model_coefs_reduced[term == "summed_tonnes_scaled_byreg",]

fishing_model_coefs_reduced[,survey_unit := grp][,summed_tonnes_scaled_byreg := condval]


fishing_model_coefs_reduced[,Slope_Direction := ifelse(summed_tonnes_scaled_byreg > 0, "Positive","Negative")]

#
fishing_model_coefs_reduced <- fishing_model_coefs_reduced[color_link, on = "survey_unit"]


#does it cross zero?
fishing_model_coefs_reduced[,significant := ifelse(lwr >0 & upr>0,T,ifelse(lwr<0 & upr<0,T,F))]

#delete all obs that are significant
fishing_model_coefs_reduced.r <- fishing_model_coefs_reduced[significant == F,]

#order table by coefficient
setorder(fishing_model_coefs_reduced, summed_tonnes_scaled_byreg)

BC_balanced_fishing_model_coefs_reduced.unique <- unique(fishing_model_coefs_reduced[,.(condval,condsd, lwr, upr, survey_unit, summed_tonnes_scaled_byreg, Slope_Direction, hex, Survey_Name_Season, significant)]) 

#extract color hexes
#year adj coef order
color_year_adj_order <-BC_balanced_fishing_model_coefs_reduced.unique[,hex]

#alphabetical order
BC_balanced_fishing_model_coefs_reduced.unique.alpha <- setorder(BC_balanced_fishing_model_coefs_reduced.unique, Survey_Name_Season)
color_alpha_order <- BC_balanced_fishing_model_coefs_reduced.unique.alpha[,hex]

ggplot() +
    geom_errorbar(data = fishing_model_coefs_reduced, aes(x = reorder(Survey_Name_Season, summed_tonnes_scaled_byreg) , y = summed_tonnes_scaled_byreg, label = Survey_Name_Season, ymin = lwr, ymax = upr), fill = "grey", width = 0) + #add confidence intervals
  geom_point(data = fishing_model_coefs_reduced, aes(x = reorder(Survey_Name_Season, summed_tonnes_scaled_byreg) , y = summed_tonnes_scaled_byreg, label = Survey_Name_Season, 
      fill = Slope_Direction), stat = 'identity', shape = 21, color = "black") +
  scale_fill_manual(values = c("white","black"), name = "Slope Direction") +
  geom_point(data = fishing_model_coefs_reduced.r, aes(x = reorder(Survey_Name_Season, summed_tonnes_scaled_byreg) , y = summed_tonnes_scaled_byreg, 
          label = Survey_Name_Season), stat = 'identity', fill = "grey",color = "grey", shape = 21) +
  geom_hline(yintercept = 0) +
  xlab("Survey unit") +
  ylab("Balanced BC dissimilarity ~\nrelative fishing pressure") +
  coord_flip() +
  theme_classic()

Warning: Ignoring unknown parameters: `fill`Warning: Ignoring unknown aesthetics: labelWarning: Ignoring unknown aesthetics: labelWarning: Ignoring unknown aesthetics: label

#scatterplot
dissimilarities_temp_fishing_regstats.na <- na.omit(dissimilarities_temp_fishing_regstats, cols = "summed_tonnes_scaled_byreg") #delete any rows with NAs

#Predicted values
dissimilarities_temp_fishing_regstats.na$pred_summed_tonnes_scaled_byreg <- predict(allreg_dissimilarity_fishing_mean_mod,re.form=NA)  ## population level
dissimilarities_temp_fishing_regstats.na$pred_summed_tonnes_scaled_byreg_individual <- predict(allreg_dissimilarity_fishing_mean_mod) ## individual level

#confidence intervals
fishing_confit <- confint(allreg_dissimilarity_fishing_mean_mod, oldNames = F)

Computing profile confidence intervals ...

#upper lower bounds from confidence intervals
dissimilarities_temp_fishing_regstats.na[,pred_summed_tonnes_scaled_byreg_lwr := fishing_confit[5,1] + fishing_confit[6,1]*summed_tonnes_scaled_byreg]#lower confidence interval
dissimilarities_temp_fishing_regstats.na[,pred_summed_tonnes_scaled_byreg_upr := fishing_confit[5,2] + fishing_confit[6,2]*summed_tonnes_scaled_byreg]#upper confidence interval


#for all regions
ggplot(data = dissimilarities_temp_fishing_regstats.na) +
  labs(x = "Relative fishing pressure",  y = "BC balanced dissimilarity") +
  geom_point(aes(x = summed_tonnes_scaled_byreg, y = bray_curtis_dissimilarity_balanced_mean, color = survey_unit), alpha = 0.3) +
  geom_line(aes(x = summed_tonnes_scaled_byreg, y = pred_summed_tonnes_scaled_byreg_individual, color = survey_unit)) +
  geom_line(aes(x = summed_tonnes_scaled_byreg, y = pred_summed_tonnes_scaled_byreg), lwd = 1) +
  scale_color_manual(values =  palette_37) +
  theme_classic()

#for all regions faceted
ggplot(data = dissimilarities_temp_fishing_regstats.na) +
  labs(x = "Relative fishing pressure",  y = "BC balanced dissimilarity") +
  geom_point(aes(x = summed_tonnes_scaled_byreg, y = bray_curtis_dissimilarity_balanced_mean, color = survey_unit), alpha = 0.3) +
  geom_line(aes(x = summed_tonnes_scaled_byreg, y = pred_summed_tonnes_scaled_byreg_individual, color = survey_unit)) +
  scale_color_manual(values =  palette_37) +
  facet_wrap(~survey_unit, scales = "free") +
  theme_classic() +
    theme(legend.position = "null")

NA
NA
NA

Other variables to consider

• Species number (spp_count_annual, different value each region and year)

Species number in a year vs. dissimilarity

ggplot(data = dissimilarities_temp_fishing_regstats) +
  labs(x = "Spp number",  y = "BC balanced dissimilarity") +
  geom_point(aes(x = spp_count_annual, y = bray_curtis_dissimilarity_balanced_mean, color = survey_unit), alpha = 0.3) +
  scale_color_manual(values =  palette_37) +
 # facet_wrap(~survey_unit, scales = "free") +
  theme_classic() +
    theme(legend.position = "null")

Species number in a year vs. dissimilarity ~ year coefficient

dissimilarities_temp_fishing_regstats.unique <- unique(dissimilarities_temp_fishing_regstats[,.(year, spp_count_annual, bray_coef, survey_unit)])

ggplot(data = dissimilarities_temp_fishing_regstats.unique) +
  labs(x = "Spp number",  y = "Bray Curtis dissimilarity ~ year coefficient") +
  geom_point(aes(x = spp_count_annual, y = bray_coef, color = survey_unit), alpha = 0.3) +
  scale_color_manual(values =  palette_37) +
  geom_hline(yintercept = 0) +
  theme_classic()

Regions that are differentiating (positive slope) tend to have fewer #s of species. The few regions with a lot of species (over 150; Gulf of Mexico, West Coast US, Northeast US, Gulf of Alaska), all have negative coefficients (homogenizing). This is the opposite of what I would have expected, as more species across a survey region would intuitively allow for more opportunities to differentiate.